{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cdc29b706b8f89a",
   "metadata": {},
   "source": [
    "# Convergence of powerflow"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1b584af942d8b5e",
   "metadata": {},
   "source": [
    "This Notebook gives an insight into how the loadflow converges and what influences it. Original analysis was done by bergkvist (https://github.com/bergkvist) in May 2020, see https://github.com/e2nIEE/pandapower/issues/786.\n",
    "\n",
    "The resolution is set low, if you are interested in higher resolution please increase the resolution, but do not commit. Since this will slow down the tests tremendously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df679f33ab734e1a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:01:54.901045Z",
     "start_time": "2025-11-25T15:01:52.740101Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandapower as pp\n",
    "import pandapower.networks as nw\n",
    "import matplotlib.pyplot as plt\n",
    "import copy\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f29f58dbd1ba1ae5",
   "metadata": {},
   "source": [
    "## Regions of convergence: scaling r and x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "initial_id",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:02:09.394577Z",
     "start_time": "2025-11-25T15:01:54.905719Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cig = nw.cigre_networks.create_cigre_network_hv()\n",
    "def convergence_test(r_ohm_factor, x_ohm_factor):\n",
    "    net = copy.deepcopy(cig)\n",
    "    net.line.r_ohm_per_km *= r_ohm_factor\n",
    "    net.line.x_ohm_per_km *= x_ohm_factor\n",
    "    try:\n",
    "        pp.runpp(net)\n",
    "        return 1\n",
    "    except Exception as e:\n",
    "        return np.nan\n",
    "\n",
    "x = np.linspace(-12, 5, 20)\n",
    "y = np.linspace(-3, 2, 20)\n",
    "z = np.array([[ convergence_test(xi, yi) for xi in x ] for yi in y ])\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "cp = ax.contourf(x, y, z)\n",
    "ax.set_xlabel('line.r_ohm_per_km *= factor')\n",
    "ax.set_ylabel('line.x_ohm_per_km *= factor')\n",
    "ax.set_title('Cigre HV: Region of convergence')\n",
    "plt.plot([0, 0], [y[0], y[-1]], '-k')\n",
    "plt.plot([x[0], x[-1]], [0, 0], '-k')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3eac2c4012d0b1cc",
   "metadata": {},
   "source": [
    "Notice the split in the middle! If x_ohm_per_km is below a threshold for this network, it doesn't matter what value we set for r_ohm_per_km, as it will always diverge."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13ca69151363758d",
   "metadata": {},
   "source": [
    "### Cigre MV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60e87d7b8a42f785",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:02:18.855865Z",
     "start_time": "2025-11-25T15:02:09.464847Z"
    }
   },
   "outputs": [],
   "source": [
    "cig = nw.cigre_networks.create_cigre_network_mv()\n",
    "def convergence_test(r_ohm_factor, x_ohm_factor):\n",
    "    net = copy.deepcopy(cig)\n",
    "    net.line.r_ohm_per_km *= r_ohm_factor\n",
    "    net.line.x_ohm_per_km *= x_ohm_factor\n",
    "    try:\n",
    "        pp.runpp(net)\n",
    "        return 1\n",
    "    except Exception as e:\n",
    "        return np.nan\n",
    "\n",
    "x = np.linspace(-200, 10, 20)\n",
    "y = np.linspace(-80, 80, 20)\n",
    "z = np.array([[ convergence_test(xi, yi) for xi in x ] for yi in y ])\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "cp = ax.contourf(x, y, z)\n",
    "ax.set_xlabel('line.r_ohm_per_km *= factor')\n",
    "ax.set_ylabel('line.x_ohm_per_km *= factor')\n",
    "ax.set_title('Cigre HV: Region of convergence')\n",
    "plt.plot([0, 0], [y[0], y[-1]], '-k')\n",
    "plt.plot([x[0], x[-1]], [0, 0], '-k')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aae363c038dc926b",
   "metadata": {},
   "source": [
    "Notice that the region of convergence is a lot larger. In this first image we can't even see the \"hole\" around the origin."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7d6595674939e8c",
   "metadata": {},
   "source": [
    "### Cigre LV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adf8bd33d9200074",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:02:28.235547Z",
     "start_time": "2025-11-25T15:02:18.874963Z"
    }
   },
   "outputs": [],
   "source": [
    "cig = nw.cigre_networks.create_cigre_network_lv()\n",
    "def convergence_test(r_ohm_factor, x_ohm_factor):\n",
    "    net = copy.deepcopy(cig)\n",
    "    net.line.r_ohm_per_km *= r_ohm_factor\n",
    "    net.line.x_ohm_per_km *= x_ohm_factor\n",
    "    try:\n",
    "        pp.runpp(net)\n",
    "        return 1\n",
    "    except Exception as e:\n",
    "        return np.nan\n",
    "\n",
    "x = np.linspace(-200, 10, 20)\n",
    "y = np.linspace(-80, 80, 20)\n",
    "z = np.array([[ convergence_test(xi, yi) for xi in x ] for yi in y ])\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "cp = ax.contourf(x, y, z)\n",
    "ax.set_xlabel('line.r_ohm_per_km *= factor')\n",
    "ax.set_ylabel('line.x_ohm_per_km *= factor')\n",
    "ax.set_title('Cigre HV: Region of convergence')\n",
    "plt.plot([0, 0], [y[0], y[-1]], '-k')\n",
    "plt.plot([x[0], x[-1]], [0, 0], '-k')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaf14b62dcd93d7b",
   "metadata": {},
   "source": [
    "### Observations\n",
    "- The powerflow calculations will converge for negative r and x-values. The region of convergence also seems to generally be bigger on the negative side for some reason. I have no intuition about what a negative resistance value means.\n",
    "- It seems like the higher the voltage, the more \"sensitive\" the network becomes to parameter scaling. (at least impedance).\n",
    "- Sometimes r or x can be 0 without causing divergence. In the case of Cigre HV, if x_ohm_per_km is below a threshold, it seems to guarantee divergence."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb4f2195b3f7e5c8",
   "metadata": {},
   "source": [
    "## vk_percent vs vkr_percent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72df4e06856ed9d6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:02:38.756882Z",
     "start_time": "2025-11-25T15:02:28.261174Z"
    }
   },
   "outputs": [],
   "source": [
    "cig = nw.cigre_networks.create_cigre_network_hv()\n",
    "def convergence_test(vk_percent, vkr_percent):\n",
    "    net = copy.deepcopy(cig)\n",
    "    net.trafo.vkr_percent = vkr_percent\n",
    "    net.trafo.vk_percent = vk_percent\n",
    "    try:\n",
    "        pp.runpp(net)\n",
    "        return 1\n",
    "    except Exception as e:\n",
    "        return np.nan\n",
    "\n",
    "x = np.linspace(0, 30, 20)\n",
    "y = np.linspace(0, 15, 20)\n",
    "z = np.array([[ convergence_test(xi, yi) for xi in x ] for yi in y ])\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "cp = ax.contourf(x, y, z)\n",
    "ax.set_xlabel('trafo.vk_percent')\n",
    "ax.set_ylabel('trafo.vkr_percent')\n",
    "ax.set_title('Cigre HV: Region of convergence')\n",
    "plt.plot([0, 0], [y[0], y[-1]], '-k')\n",
    "plt.plot([x[0], x[-1]], [0, 0], '-k')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "268e4fe1d17a3ade",
   "metadata": {},
   "source": [
    "Seems a bit more well behaved than the impedance-ROCs. If this is actually a perfect triangle, then maybe we could predict based on vkr_percent and vk_percent values whether the powerflow will diverge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a92c271008f9ecdf",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:03:01.654936Z",
     "start_time": "2025-11-25T15:02:38.782598Z"
    }
   },
   "outputs": [],
   "source": [
    "cig = nw.cigre_networks.create_cigre_network_hv()\n",
    "def convergence_test(vk_percent, vkr_percent):\n",
    "    net = copy.deepcopy(cig)\n",
    "    net.trafo.vkr_percent = vkr_percent\n",
    "    net.trafo.vk_percent = vk_percent\n",
    "    try:\n",
    "        pp.runpp(net)\n",
    "        return 1\n",
    "    except Exception as e:\n",
    "        return np.nan\n",
    "\n",
    "x = np.linspace(-50, 50, 30)\n",
    "y = np.linspace(-50, 12, 30)\n",
    "z = np.array([[ convergence_test(xi, yi) for xi in x ] for yi in y ])\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "cp = ax.contourf(x, y, z)\n",
    "ax.set_xlabel('trafo.vk_percent')\n",
    "ax.set_ylabel('trafo.vkr_percent')\n",
    "ax.set_title('Cigre HV: Region of convergence')\n",
    "plt.plot([0, 0], [y[0], y[-1]], '-k')\n",
    "plt.plot([x[0], x[-1]], [0, 0], '-k')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a70d273c3740e66",
   "metadata": {},
   "source": [
    "### Exploring vk_percent and vkr_percent convergence a bit more"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d382965a157552c7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:04:08.294409Z",
     "start_time": "2025-11-25T15:03:01.676848Z"
    }
   },
   "outputs": [],
   "source": [
    "cig = nw.cigre_networks.create_cigre_network_mv()\n",
    "\n",
    "def convergence_test(vk_percent, vkr_percent):\n",
    "    net = copy.deepcopy(cig) # Deepcopy is much faster than loading the network again\n",
    "    net.trafo.vk_percent = vk_percent\n",
    "    net.trafo.vkr_percent = vkr_percent\n",
    "    try:\n",
    "        pp.runpp(net)\n",
    "        return 1\n",
    "    except Exception as e:\n",
    "        return np.nan\n",
    "\n",
    "# Binary search to find the limit between convergence and divergence for vkr_percent\n",
    "def find_vkr_percent_limit(vk_percent, vkr_conv, vkr_div, tolerance=1e-3):\n",
    "    if abs(vkr_div - vkr_conv) < tolerance: return vkr_conv\n",
    "    vkr_test = 0.5 * (vkr_div + vkr_conv)\n",
    "    if np.isnan(convergence_test(vk_percent, vkr_test)):\n",
    "        return find_vkr_percent_limit(vk_percent, vkr_conv, vkr_test, tolerance)\n",
    "    else:\n",
    "        return find_vkr_percent_limit(vk_percent, vkr_test, vkr_div, tolerance)\n",
    "\n",
    "# Assuming vkr_percent=0 converges, vk_percent=40 diverges, the limits will be found:\n",
    "vk_percent = np.linspace(0, 50, 151)\n",
    "vkr_percent = np.vectorize(find_vkr_percent_limit)(vk_percent, vkr_conv=0.0, vkr_div=40.0)\n",
    "plt.plot(vk_percent, vkr_percent)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34ffc88ee7fbe8c3",
   "metadata": {},
   "source": [
    "Notice that up to some value for vk_percent, 0 <= vkr_percent <= vk_percent will cause converge. But we do actually get divergence here when vkr_percent > vk_percent.\n",
    "\n",
    "The trailing 0s in the end of the plot means that no solution between 0 and 40 was found."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f617f5b36bd73798",
   "metadata": {},
   "source": [
    "## As for c_nf_per_km, and the effect of scaling the voltage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0fd641321360e40",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:04:35.726863Z",
     "start_time": "2025-11-25T15:04:08.345859Z"
    }
   },
   "outputs": [],
   "source": [
    "cig = nw.cigre_networks.create_cigre_network_hv()\n",
    "def convergence_test(vn_kv_factor, c_nf_factor):\n",
    "    net = copy.deepcopy(cig)\n",
    "    net.bus.vn_kv *= vn_kv_factor\n",
    "    net.trafo.vn_hv_kv *= vn_kv_factor\n",
    "    net.trafo.vn_lv_kv *= vn_kv_factor\n",
    "    net.shunt.vn_kv *= vn_kv_factor\n",
    "    net.line.c_nf_per_km *= c_nf_factor\n",
    "    try:\n",
    "        pp.runpp(net)\n",
    "        return 1\n",
    "    except Exception as e:\n",
    "        return np.nan\n",
    "\n",
    "x = np.linspace(0, 20, 30)\n",
    "y = np.linspace(-10, 10, 30)\n",
    "z = np.array([[ convergence_test(xi, yi) for xi in x ] for yi in y ])\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "cp = ax.contourf(x, y, z)\n",
    "ax.set_xlabel('factor: bus.vn_kv')\n",
    "ax.set_ylabel('factor: line.c_nf_per_km')\n",
    "ax.set_title('Cigre HV: Region of convergence')\n",
    "plt.plot([0, 0], [y[0], y[-1]], '-k')\n",
    "plt.plot([x[0], x[-1]], [0, 0], '-k')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72a86a81eceeb4d5",
   "metadata": {},
   "source": [
    "These regions of convergence look surprisingly weird/interesting."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a2082a7c4336bd5",
   "metadata": {},
   "source": [
    "## Some more analysis on r_ohm and x_ohm-scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb609fb04096560d",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:05:20.899102Z",
     "start_time": "2025-11-25T15:04:35.748148Z"
    }
   },
   "outputs": [],
   "source": [
    "def from_polar(r, th):\n",
    "    return r * np.cos(th), r * np.sin(th)\n",
    "\n",
    "@np.vectorize\n",
    "def limit_polar_log_search(test_fn, theta, log_radius_conv, log_radius_div,\n",
    "                           tolerance=1e-3, max_iterations=50):\n",
    "    # Check that radius_conv does not cause divergence\n",
    "    if np.isnan(test_fn(*from_polar(np.exp(log_radius_conv), theta))):\n",
    "        return np.nan\n",
    "\n",
    "    # Check that radius_div does not cause convergence\n",
    "    if test_fn(*from_polar(np.exp(log_radius_div), theta)) == 1:\n",
    "        return np.nan\n",
    "\n",
    "    while (np.abs(log_radius_div - log_radius_conv) > tolerance\n",
    "           and max_iterations > 0):\n",
    "        max_iterations -= 1\n",
    "        log_radius_test = (log_radius_div + log_radius_conv) / 2\n",
    "        if np.isnan(test_fn(*from_polar(np.exp(log_radius_test), theta))):\n",
    "            log_radius_div = log_radius_test\n",
    "        else:\n",
    "            log_radius_conv = log_radius_test\n",
    "\n",
    "    return np.exp(log_radius_conv)\n",
    "\n",
    "cig = nw.cigre_networks.create_cigre_network_hv()\n",
    "def convergence_test(r_ohm_factor, x_ohm_factor):\n",
    "    net = copy.deepcopy(cig) # Deepcopy is much faster than loading the network again\n",
    "    net.line.r_ohm_per_km *= r_ohm_factor\n",
    "    net.line.x_ohm_per_km *= x_ohm_factor\n",
    "    try:\n",
    "        pp.runpp(net)\n",
    "        return 1\n",
    "    except Exception as e:\n",
    "        return np.nan\n",
    "\n",
    "theta = np.linspace(0, np.pi/2, 50)\n",
    "radius_upper = limit_polar_log_search(convergence_test, theta, log_radius_conv=0, log_radius_div=50)\n",
    "radius_lower = limit_polar_log_search(convergence_test, theta, log_radius_conv=0, log_radius_div=-50)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(theta, radius_lower)\n",
    "ax.plot(theta, radius_upper)\n",
    "ax.set_xlabel('impedance scaling factor: theta')\n",
    "ax.set_ylabel('impedance scaling factor: radius')\n",
    "ax.set_title('Cigre HV: Upper and lower bounds for convergence')\n",
    "ax.set_yscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1737077de342d8cc",
   "metadata": {},
   "source": [
    "### Observations\n",
    "- Now we are able to see the upper and lower limits in a single visualization. This is able to tell us something about how robust our values for r_ohm and x_ohm are. Essentially: how safe is it to scale them up or down.\n",
    "- As we saw in the previous 2D-ROC-plot, a small angle will guarantee divergence for Cigre HV. This is kind of hard to notice in the plot (where nan-values are not shown). This should probably be highlighted in some kind of way.\n",
    "- This is a lot more efficient than the 2D-ROC-pixel-plots for finding upper and lower bounds with high accuracy.\n",
    "- The slope of the upper bound tells us whether the network is more resilient to r_ohm increasing or x_ohm increasing. Maybe this can tell us something interesting about the network.\n",
    "- The lower bound is very \"uneven\" while the upper bound seems smooth. The upper and lower bounds seem to have a slope in the same direction."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b12ccd8a7d7fa729",
   "metadata": {},
   "source": [
    "## Voltage stability\n",
    "\n",
    "ext_grid(1pu) <--> bus(10kV) <--> line <--> bus(10kV) <--> load(1MW)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5101302641930470",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:05:23.065018Z",
     "start_time": "2025-11-25T15:05:20.934593Z"
    }
   },
   "outputs": [],
   "source": [
    "net = pp.create_empty_network()\n",
    "b1 = pp.create_bus(net, vn_kv=10.0, name='b1')\n",
    "b2 = pp.create_bus(net, vn_kv=10.0, name='b2')\n",
    "pp.create_load(net, bus=b2, p_mw=1.)\n",
    "pp.create_ext_grid(net, bus=b1, vm_pu=1.0)\n",
    "pp.create_line_from_parameters(net, from_bus=b1, to_bus=b2, length_km=1,\n",
    "                                   r_ohm_per_km=1, x_ohm_per_km=0.1,\n",
    "                                   c_nf_per_km=0.1, max_i_ka=10)\n",
    "\n",
    "def convergence_test(c_nf_factor):\n",
    "    nt = copy.deepcopy(net)\n",
    "    nt.line.c_nf_per_km *= c_nf_factor\n",
    "    try:\n",
    "        pp.runpp(nt)\n",
    "        return nt.res_bus.vm_pu\n",
    "    except Exception as e:\n",
    "        return np.array([0., 0.])\n",
    "\n",
    "x = np.linspace(1e4, 3e8, 100)\n",
    "z = np.array([ convergence_test(xi) for xi in x ])\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.xlabel('line.c_nf_per_km')\n",
    "plt.ylabel('bus.vm_pu')\n",
    "fig.suptitle('ext_grid(1pu) <--> bus(10kV) <--> line <--> bus(10kV) <--> load(1MW)')\n",
    "plt.plot(x, z, 'o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5fa31e164135575",
   "metadata": {},
   "source": [
    "### Observations\n",
    "- The blue line represents the bus connected to the slack bus, so it makes sense that the voltage here remains constant. Notice the \"holes\" in this line, however - showing where the powerflow diverges.\n",
    "- The orange line represents the load bus, at the other end of the line. The voltage starts rising, before diverging randomly and fluctuating between ~0.2 and 1.0-1.05.\n",
    "- Since Netwton-Rhapson can only converge to a single solution, maybe this means we actually have 2 solutions here, and it is unpredictable which one we will converge to? And one could guess this region might also be highly unlinear, meaning the solution can easily diverge by hitting a bump on the way to a solution.\n",
    "- Based on the 2D-ROC plots we have seen that the ROC sometimes seem to become fractal-like at the border (although in some places it is actually very smooth).\n",
    "\n",
    "### Questions\n",
    "- Does the situation with vm_pu = 0.2 in the result correspond to a \"voltage collapse\", like those that have caused several major power grid blackouts?\n",
    "- When there are multiple solutions, how would reality \"pick one of them\"?\n",
    "    My guess is that both would be valid steady-states of the grid in reality. But in this region the network would be extremely sensitive to disturbances - that could throw it into a voltage-collapsed state.\n",
    "- Why does it seem like there is a continuous transition to the collapsed state?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "926ebfafaa91022",
   "metadata": {},
   "source": [
    "## Divergence doesn't neccesarily correspond to voltage collapse!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2641d5a3a50ffe6c",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:05:35.599961Z",
     "start_time": "2025-11-25T15:05:23.091451Z"
    }
   },
   "outputs": [],
   "source": [
    "def binary_search(fn, x_ok, x_err, tolerance=1e-6, max_iterations=20):\n",
    "    def try_fn(x):\n",
    "        try:\n",
    "            fn(x)\n",
    "            return 1\n",
    "        except:\n",
    "            return 0\n",
    "\n",
    "    assert try_fn(x_ok) == 1\n",
    "    assert try_fn(x_err) == 0\n",
    "\n",
    "    while abs(x_err - x_ok) > tolerance and max_iterations > 0:\n",
    "        x_guess = (x_ok + x_err) / 2\n",
    "        if try_fn(x_guess) == 0:\n",
    "            x_err = x_guess\n",
    "        else:\n",
    "            x_ok = x_guess\n",
    "        max_iterations -= 1\n",
    "\n",
    "    return x_ok\n",
    "\n",
    "def create_network(p_mw):\n",
    "    net = pp.create_empty_network()\n",
    "    b0, b1 = pp.create_buses(net, nr_buses=2, vn_kv=10)\n",
    "    pp.create_ext_grid(net, bus=b0, vm_pu=1)\n",
    "    pp.create_line_from_parameters(net, from_bus=b0, to_bus=b1, length_km=1,\n",
    "                                   r_ohm_per_km=1, x_ohm_per_km=0.1,\n",
    "                                   c_nf_per_km=0, max_i_ka=10)\n",
    "    pp.create_load(net, bus=b1, p_mw=p_mw)\n",
    "    return net\n",
    "\n",
    "@np.vectorize\n",
    "def find_voltage(p_mw, vm_pu_init):\n",
    "    net = create_network(p_mw)\n",
    "    try:\n",
    "        pp.runpp(net, init_vm_pu=[1, vm_pu_init])\n",
    "        return net.res_bus.iloc[1].vm_pu\n",
    "    except:\n",
    "        return np.nan\n",
    "\n",
    "# Find the maximum value for p_mw where the network converges\n",
    "p_mw_max = binary_search(lambda p_mw: pp.runpp(create_network(p_mw), init='flat'), 0, 50)\n",
    "\n",
    "# Show that in there are two different solutions for every converging p_mw\n",
    "p_mw = np.linspace(0, p_mw_max, 50)\n",
    "vm_pu_high = find_voltage(p_mw, vm_pu_init=1)\n",
    "vm_pu_low = find_voltage(p_mw, vm_pu_init=0.1)\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.suptitle('Bifurcation diagram (r_ohm=1, x_ohm=0.1, c_nf=0)')\n",
    "plt.plot(p_mw, vm_pu_high)\n",
    "plt.plot(p_mw, vm_pu_low)\n",
    "plt.xlabel('p_mw')\n",
    "plt.ylabel('vm_pu')\n",
    "plt.legend(['default solution','alternative solution'])\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcb53b6621627f9b",
   "metadata": {},
   "source": [
    "Divergence is caused by trying to use more power than the maximum power transfer theorem allows.\n",
    "https://en.wikipedia.org/wiki/Maximum_power_transfer_theorem\n",
    "\n",
    "The two solutions simply correspond to the two possible load impedances that yield the same power consumption."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44ac0d7009da769f",
   "metadata": {},
   "source": [
    "## Setting c_nf_per_km to a high value (1e7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ffcc01cb5d9012b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:05:48.316688Z",
     "start_time": "2025-11-25T15:05:35.626481Z"
    }
   },
   "outputs": [],
   "source": [
    "def binary_search(fn, x_ok, x_err, tolerance=1e-6, max_iterations=20):\n",
    "    def try_fn(x):\n",
    "        try:\n",
    "            fn(x)\n",
    "            return 1\n",
    "        except:\n",
    "            return 0\n",
    "\n",
    "    assert try_fn(x_ok) == 1\n",
    "    assert try_fn(x_err) == 0\n",
    "\n",
    "    while abs(x_err - x_ok) > tolerance and max_iterations > 0:\n",
    "        x_guess = (x_ok + x_err) / 2\n",
    "        if try_fn(x_guess) == 0:\n",
    "            x_err = x_guess\n",
    "        else:\n",
    "            x_ok = x_guess\n",
    "        max_iterations -= 1\n",
    "\n",
    "    return x_ok\n",
    "\n",
    "def create_network(p_mw):\n",
    "    net = pp.create_empty_network()\n",
    "    b0, b1 = pp.create_buses(net, nr_buses=2, vn_kv=10)\n",
    "    pp.create_ext_grid(net, bus=b0, vm_pu=1)\n",
    "    pp.create_line_from_parameters(net, from_bus=b0, to_bus=b1, length_km=1,\n",
    "                                   r_ohm_per_km=1, x_ohm_per_km=0.1,\n",
    "                                   c_nf_per_km=1e7, max_i_ka=10)\n",
    "    pp.create_load(net, bus=b1, p_mw=p_mw)\n",
    "    return net\n",
    "\n",
    "@np.vectorize\n",
    "def find_voltage(p_mw, vm_pu_init):\n",
    "    net = create_network(p_mw)\n",
    "    try:\n",
    "        pp.runpp(net, init_vm_pu=[1, vm_pu_init])\n",
    "        return net.res_bus.iloc[1].vm_pu\n",
    "    except:\n",
    "        return np.nan\n",
    "\n",
    "# Find the maximum value for p_mw where the network converges\n",
    "p_mw_max = binary_search(lambda p_mw: pp.runpp(create_network(p_mw), init='flat'), 0, 50)\n",
    "\n",
    "# Show that in there are two different solutions for every converging p_mw\n",
    "p_mw = np.linspace(0, p_mw_max, 50)\n",
    "vm_pu_high = find_voltage(p_mw, vm_pu_init=1)\n",
    "vm_pu_low = find_voltage(p_mw, vm_pu_init=0.1)\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.suptitle('Bifurcation diagram (r_ohm=1, x_ohm=0.1, c_nf=1e7)')\n",
    "plt.plot(p_mw, vm_pu_high)\n",
    "plt.plot(p_mw, vm_pu_low)\n",
    "plt.xlabel('p_mw')\n",
    "plt.ylabel('vm_pu')\n",
    "plt.legend(['default solution','alternative solution'])\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9b5d43da61410cd",
   "metadata": {},
   "source": [
    "Some things to notice here:\n",
    "\n",
    "- Maximum power transfer at vm_pu ~0.3, p_mw ~17.5 MW\n",
    "- The diagram is not symmetric. The lower part is now taller than the upper part.\n",
    "- Low voltage (alternative) solutions suffer from numerical instability. (the alternative solution corresponds to the load impedance being close to 0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fdfdec537376f2",
   "metadata": {},
   "source": [
    "## Setting x_ohm = 3 * r_ohm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a2db626484435d6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:06:01.989506Z",
     "start_time": "2025-11-25T15:05:48.344601Z"
    }
   },
   "outputs": [],
   "source": [
    "def binary_search(fn, x_ok, x_err, tolerance=1e-6, max_iterations=20):\n",
    "    def try_fn(x):\n",
    "        try:\n",
    "            fn(x)\n",
    "            return 1\n",
    "        except:\n",
    "            return 0\n",
    "\n",
    "    assert try_fn(x_ok) == 1\n",
    "    assert try_fn(x_err) == 0\n",
    "\n",
    "    while abs(x_err - x_ok) > tolerance and max_iterations > 0:\n",
    "        x_guess = (x_ok + x_err) / 2\n",
    "        if try_fn(x_guess) == 0:\n",
    "            x_err = x_guess\n",
    "        else:\n",
    "            x_ok = x_guess\n",
    "        max_iterations -= 1\n",
    "\n",
    "    return x_ok\n",
    "\n",
    "def create_network(p_mw):\n",
    "    net = pp.create_empty_network()\n",
    "    b0, b1 = pp.create_buses(net, nr_buses=2, vn_kv=10)\n",
    "    pp.create_ext_grid(net, bus=b0, vm_pu=1)\n",
    "    pp.create_line_from_parameters(net, from_bus=b0, to_bus=b1, length_km=1,\n",
    "                                   r_ohm_per_km=1, x_ohm_per_km=3,\n",
    "                                   c_nf_per_km=0, max_i_ka=10)\n",
    "    pp.create_load(net, bus=b1, p_mw=p_mw)\n",
    "    return net\n",
    "\n",
    "@np.vectorize\n",
    "def find_voltage(p_mw, vm_pu_init):\n",
    "    net = create_network(p_mw)\n",
    "    try:\n",
    "        pp.runpp(net, init_vm_pu=[1, vm_pu_init])\n",
    "        return net.res_bus.iloc[1].vm_pu\n",
    "    except:\n",
    "        return np.nan\n",
    "\n",
    "# Find the maximum value for p_mw where the network converges\n",
    "p_mw_max = binary_search(lambda p_mw: pp.runpp(create_network(p_mw), init='flat'), 0, 50)\n",
    "\n",
    "# Show that in there are two different solutions for every converging p_mw\n",
    "p_mw = np.linspace(0, p_mw_max, 50)\n",
    "vm_pu_high = find_voltage(p_mw, vm_pu_init=1)\n",
    "vm_pu_low = find_voltage(p_mw, vm_pu_init=0.1)\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.suptitle('Bifurcation diagram (r_ohm=1, x_ohm=3, c_nf=0)')\n",
    "plt.plot(p_mw, vm_pu_high)\n",
    "plt.plot(p_mw, vm_pu_low)\n",
    "plt.xlabel('p_mw')\n",
    "plt.ylabel('vm_pu')\n",
    "plt.legend(['default solution','alternative solution'])\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "292aece4ce1a0e6a",
   "metadata": {},
   "source": [
    "- Maximum power transfer at vm_pu ~0.6, p_mw ~12MW\n",
    "- The diagram is not symmetric.\n",
    "- Low voltage (alternative) solutions suffer from numerical instability."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc0205bccc8dab37",
   "metadata": {},
   "source": [
    "## Purely reactive line impedance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "807ab9706b6bb718",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-11-25T15:06:15.535578Z",
     "start_time": "2025-11-25T15:06:02.017914Z"
    }
   },
   "outputs": [],
   "source": [
    "def binary_search(fn, x_ok, x_err, tolerance=1e-6, max_iterations=20):\n",
    "    def try_fn(x):\n",
    "        try:\n",
    "            fn(x)\n",
    "            return 1\n",
    "        except:\n",
    "            return 0\n",
    "\n",
    "    assert try_fn(x_ok) == 1\n",
    "    assert try_fn(x_err) == 0\n",
    "\n",
    "    while abs(x_err - x_ok) > tolerance and max_iterations > 0:\n",
    "        x_guess = (x_ok + x_err) / 2\n",
    "        if try_fn(x_guess) == 0:\n",
    "            x_err = x_guess\n",
    "        else:\n",
    "            x_ok = x_guess\n",
    "        max_iterations -= 1\n",
    "\n",
    "    return x_ok\n",
    "\n",
    "def create_network(p_mw):\n",
    "    net = pp.create_empty_network()\n",
    "    b0, b1 = pp.create_buses(net, nr_buses=2, vn_kv=10)\n",
    "    pp.create_ext_grid(net, bus=b0, vm_pu=1)\n",
    "    pp.create_line_from_parameters(net, from_bus=b0, to_bus=b1, length_km=1,\n",
    "                                   r_ohm_per_km=0, x_ohm_per_km=1,\n",
    "                                   c_nf_per_km=0, max_i_ka=10)\n",
    "    pp.create_load(net, bus=b1, p_mw=p_mw)\n",
    "    return net\n",
    "\n",
    "@np.vectorize\n",
    "def find_voltage(p_mw, vm_pu_init):\n",
    "    net = create_network(p_mw)\n",
    "    try:\n",
    "        pp.runpp(net, init_vm_pu=[1, vm_pu_init])\n",
    "        return net.res_bus.iloc[1].vm_pu\n",
    "    except:\n",
    "        return np.nan\n",
    "\n",
    "# Find the maximum value for p_mw where the network converges\n",
    "p_mw_max = binary_search(lambda p_mw: pp.runpp(create_network(p_mw), init='flat'), 0, 50)\n",
    "\n",
    "# Show that in there are two different solutions for every converging p_mw\n",
    "p_mw = np.linspace(0, p_mw_max, 50)\n",
    "vm_pu_high = find_voltage(p_mw, vm_pu_init=1)\n",
    "vm_pu_low = find_voltage(p_mw, vm_pu_init=0.1)\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.suptitle('Bifurcation diagram (r_ohm=0, x_ohm=1, c_nf=0)')\n",
    "plt.plot(p_mw, vm_pu_high)\n",
    "plt.plot(p_mw, vm_pu_low)\n",
    "plt.xlabel('p_mw')\n",
    "plt.ylabel('vm_pu')\n",
    "plt.legend(['default solution','alternative solution'])\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d25bd3af2abcae10",
   "metadata": {},
   "source": [
    "- Maximum power transfer at vm_pu ~0.7, p_mw ~50MW\n",
    "- The diagram is not symmetric.\n",
    "- Low voltage (alternative) solutions suffer from numerical instability."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
